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ABSTRACT 

An investigation of a method for finding an optimum com- 
pensator for use in an autopilot is described. Given a struc- 
ture with a number of real poles and zeros v;hich is steady- 
state decoupled, the pole-zero locations are chosen to give 
a system response as close as possible to a desired response. 
The desired response can be any one or a combination of sys- 
tem outputs. The computer is used to optimize the pole-zero 
locations by means of a function minimization program. 
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INTRODUCTION 



The autopilot is used for aircraft control to hold some 
variable of the aircraft's motion constant, such as speed, 
altitude or glide slope angle, in the presence of distur- 
bances. These disturbances can originate outside the air- 
craft, e.g. random wind gusts, or within the system such as 
electrical noise in actuators and control system components. 

The problem becomes more complex when the aircraft is 
viewed as a multivariable system, i.e. more than one input 
and output. For example, a movement of the elevator control 
will change the pitch angle and the speed, and a change in 
thrust likewise. It would be desirable for precise control 
to have one input cause a change in a corresponding output 
and in no other. This condition would require that the 
transfer function matrix of the entire system be a diagonal 
matrix, which allows much less freedom in designing for sta- 
bility . 

This condition can be improved by requiring that the 
non-change in other variables (decoupling) occur only after 
the initial transients have died out, that is, once steady 
state has been reached. 

Methods have been devised to obtain steady state de- 
coupling by linear state-variable feedback and by cascade 
compensation with unity feedback (Huang) . 

In particular, for the cascade compensation case, which 
is probably more practical than state variable feedback., Huang 
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found that the rank condition for the plant matrix is no 
longer necessary as it was in the state variable feedback 
case and that pure integrators in the compensators and poles 
at the origin in the plant transfer function were helpful 
in providing steady-state decoupling. 

However, when it comes to designing the remainder of 
the compensators for stability or a desired transient res- 
ponse, one is left v;ith a tedious root-locus design process 
which to some extent requires considerable experience in 
this technique. The problem is aggravated by the fact that 
the root loci for each channel are interdependent and this 
complexity increases with the number of inputs and outputs. 
Even when completed, there is no guarantee that the design 
is an optimum one. 

The intent of this thesis is to develop a method of ap- 
plying the power of the modern digital computer and the tech- 
niques of linear programming to the process of compensator 
design. Essentially the process consists of simulating the 
system in the time domain with the poles and zeros of the 
compensators as unknowns. By having a function minimization 
program minimize the integral-squared-difference between the 
response of the simulated system and a desired response, the 
required poles and zeros are determined. 

Cantalapiedra used a similar method with the same func- 
tion minimization program to determine the optimum number of 
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poles and zeros to be used in a model of a high order system. 
In hi.i' s thesis, the integral squared difference between the 
system and model response curves was minimized for various 
configurations of the model and this so-called cost function 
served as a measure of how well the low-order model approx- 
imated the high-order system. 

Lima, in his thesis, used the same function minimization 
program to obtain the gains for a compensated ship controller 
in a model of two ships in underway replenishment. In this 
case a desired track was to be followed by one ship with 
.reference to the replenishing ship in the presence of surge 
and sway forces affecting both ships. 
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II. DESCRIPTION OF THE SYSTEM 



In this thesis, the longitudinal motion of the C-8A 
Augmentor Wing "Buffalo” aircraft presently under test and 
development at the NASA TVmes Research Center was modelled. 
There are two reasons for this choice. First, a STOL air- 
craft such as the augmentor wing Buffalo, in the landing 
or take-off phase, operates at low speeds and steep flight 
path angles and in such flight regimes a strong artificial 
stability must be applied just to make the aircraft flyable. 
Second, a model of this same aircraft was used by Huang in 
his investigation of decoupling referenced earlier; thus the 
plant equations were readily available as well as the com- 
pensator transfer functions which he determined by the root- 
locus method. 

Only the longitudinal motion of the aircraft was in- 
vestigated since the principles arrived at for this mode 
would be equally applicable to the lateral-directional mode 
in a more complicated model. The aircraft is initially con- 
sidered as flying straight and level at a speed Vq = 126.7 
ft/sec when it is subjected to a reference pitch angle step 
input of -0.25 radians at time = 0.0. These conditions 
apply to all aspects of the investigation which follows. 
Ideally, if the system were perfectly decoupled, the pitch 
angle would change to -0.25 radians and the speed and angle 
of attack would remain at their trim conditions. However, 
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as will be noted, there is a short period of transients 
which follows the initial step input which die out leaving 
speed and angle of attack at their trim conditions and the 
pitch angle at -0,25 radians. 

Figure II-l shows a schematic diagram of the system model 
consisting of the plant dynamics and the actuator dynamics and 
cascade compensators; one for each channel. The three inputs 
to the plant are 6^ , % thrust (trim = 100%); 6^ , flap an- 

gle deflection from trim in radians and '5^ , elevator deflec- 
tion from trim in radians. The three outputs from the plant 
are u , the speed variation from trim in ft/sec, a , the an- 
gle of attack in radians, and 9, the pitch angle in radians. 

The feedback of outputs to inputs is to a certain extent ar- 
bitrary, however by inspection one could feedback an output 
to the input which had the greatest effect on it. In any case, 
the scheme adopted here is that used by Huang in his thesis. 

Figure II-2 shows a schematic diagram of the system 
model as designed by Huang. The compensator poles and zeros 
which he determined by root-locus techniques, the actuator 
transfer functions and the plant equations are also shown. 

The compensator poles at the origin, necessary for decoupling, 
should be noted. The response of this model to a pitch angle 
step input of -0.25 radians is shown in Figures II-3, II-4 
and II-5 for speed, angle of attack and pitch angle respec- 
tively. These responses are used as reference responses 
during the course of the investigation. It will be noted 
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that these responses are somewhat oscillatory and of large 
amplitude. One major intent of this investigation is to 
improve (i.e., reduce) the frequency, magnitude and duration 
of these transient oscillations. 

It is assumed for this investigation that the structure 
of the compensators is fixed; that is, the numbers of poles 
and zeros in each compensator will be as designed by Huang. 
This is not to say that this is an optimum number. Indeed 
a greater number of poles and zeros in each compensator may 
improve the response over what is possible with the present 
configuration. A lesser degree of complexity and an inter- 
est in determining how much could be done with the present 
configuration dictated the choice. Similarly, the gains in 
each compensator are considered as fixed at the values shovm 
in Figure II-2, though there is no reason why they could not 
be made a. further three variables. 
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III. TECHNIQUE FOR SOLUTION 



A. THE COST FUNCTION 

The point in question for the investigation is to deter- 
mine whether and at what cost a given 3X3 plant can be made 
to follow a desired output response by determining the poles 
and zeros of cascade compensators of a given configuration 
using the digital computer. When given a curve of the de- 
sired output response, the most obvious measure of hov; 
closely the model has matched it, is the integral-squared- 
difference between the desired output and the actual output 
between t = 0.0 and some final time t^ . This so-called 
cost function J has the form: 

^f ^ 

^ w. (y , . - y . ) ^dt 
^ X -^di -^ai 

i=l 

where y^^ = the desired i th output at time t 
= the actual i th output at time t 
w^ = weighting factor for i th output 

The w^ can be chosen to penalize a more important out- 
put deviation over a less important one or, as in the case 
of the present investigation, to give equal weight to all 
outputs which vary due to differences in units. It is not 
necessary that the cost function be a function of all outputs. 
It could be a function of the one of immediate importance 
with subsequent checking of the other outputs to follow. 
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Both types of cost functions were used in this investigation. 

The choice of t^ is an important one. If t^ is too 
long, the cost of the calculations in computer time becomes 
prohibitive, especially if the time steps used in integration 
are small. On the other hand if t^ is too short, the re- 
sults may indicate a good following of the desired response 
for t < t^, but the system may diverge and go unstable for 
t > t^. A compromise was reached for this investigation 
choosing tv/o times; a t^ = 25 seconds, the approximate 
settling time of the original theta response (see Figure II-5) , 
and t^ =12.5 seconds. The calculations for curve matching 
were done over the shorter period and the results checked by 
simulating the system and plotting the responses over the 
longer period. 

Thus it will be seen that J is a direct measure of how 
well the sytem output response curve matches a desired output 
response curve over the simulation time chosen. If J reaches 
a minimum it can be assumed that, barring the chance of a 
local rather than a global minimum, the "best" combination 
of poles and zeros has been found. Alternately, even though 
a minimum has not been found, a visual comparison of the two 
responses may indicate that the poles and zeros determined 
are close enough to an optimum set. 

B. THE FUNCTION MINIMIZATION SUBROUTINE 

1 . General 

The function minimization subroutine chosen for this 
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investigation is one which has seen some use at the Naval 
Postgraduate School in a variety of disciplines. At the 
time the present work was proceeding it was undergoing the 
incorporation of some improvements and the author had the 
benefit of this more up-to-date version. 

The subroutine is called BOXPLX and was programmed 
by R. R. Hilleary of the Naval Postgraudate School Computer 
Center. It will find the minimum of any arbitrary function, 
linear or non-linear subject to explicit constraints on the 
variables or implicit constraints on functions of the vari- 
ables. It will handle a maximum of 25 variables., independent 
and auxiliary. With the subroutine must be appended two 
user-generated function subroutines to calculate the cost 
function and check for implicit constraint violations. A 
brief description of BOXPLX is included at Appendix A to 
allow a better understanding of some of the processes which 
follow and to guide further workers in this area. 

2 . Constraints and Start Points 

The variables in BOXPLX are allowed to move within 
a feasible region (n-dimensional space, where n = number of 
variables) defined by upper and lower bounds on their values. 
These values are given as input to BOXPLX for any given pro- 
blem. 

It was assumed that all poles and zeros would be in 
the left half plane for reasons of stability and thus an 
absolute lower bound of 0.0 was set, but it may be necessary 
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or desireable to set some positive value as a lower bound 
depending on the results of a given run. 

The choice of upper bounds is more difficult. If 
they are set too high the feasible region in which BOXPLX 
searches for a minimum becomes large with a consequent 
heavy usage of computer time in reaching a minimum. If 
the bounds are too low, the "best" answer for a given run 
may be right on the bounds, necessitating a relaxation and 
further runs to proceed toward the minimum cost function. 

Some judgment is necessary in choosing these upper and lov/er 
bounds . 

The choice of the starting values for the variables 
is no less critical in avoiding wastage of computer time. 

A preliminary root-locus investigation of the system may 
give some idea of starting values which will be close to the 
"best" values, but it was decided that the real power of' 
the method described herein v/ould be in its application to 
the case where one had no idea where the ansv/ers were. Thus, 
the criteria mentioned above were not brought to bear too 
heavily on the problem. 

3 . The Function and Constraint Evaluation Sub-Functions 

The function evaluation sub-function, FE(X,XDATA), 

is the routine which actually calculates the cost-function. 

For every set (vertex) of poles and zeros generated by BOXPLX 
% 

this routine simulates the system response comparing it at 
every time step with the desired response, subtracting the 
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difference, squaring it and integrating the result from 
t = 0.0 to t = t^ . 

KE(X) checks to see if there are any implicit con- 
straint violations in the generated vertex. Implicit con- 
straints may be any arbitrary function of the variables the 
user may desire, e.g. the product of tv;o poles must be less 
than a given number. In this investigation no implicit con- 
straints were used. 

4 . The Desired Response (XDATA) 

The response curve against which the system, with its 
current set of poles and zeros, is compared may be read into 
the programme as a set of data cards defining the amplitude 
of any arbitrary response at each time step of the simulation 
interval. Thus a curve which lias no deterministic equation 
may be used to define a desired response. 

If the equation of a desired response curve is known, 
this equation could be included in function FE and the ampli- 
tudes computed simultaneously with the actual system response. 

In this investigation, data cards of the desired re- 
sponse were used, as being the most general method of spe- 
cifying the response. 

5 . Subroutine Graph 

This subroutine is appended to the main programme to 
plot the desired and actual response on a single graph when 
exit occurs from BOXPLX with the "best" set of poles and 
zeros. In. addition, it can be made to plot other responses 
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which are not being compared with a desired response, to 
check on their quality. 

C. CONVERSION OF SYSTEM EQUATIONS TO STATE VARIABLE FORM 

1 . General 

The system had previously been modeled by Huang 
using DSL (Digital Simulation Language) for which a program 
deck was available. However, since DSL is a language on its 
own, it was not possible to use it in its existing form; in 
an interactive mode in a Fortran program with the function 
minimization sub-routine. Thus, it was decided to first 
convert the system equations to state variable form so that 
they could be integrated by a built-in integration routine. 

Since the poles and zeros of the compensators will 
be the unknown parameters and it is intended to use the same 
quantity and configuration of these parameters as in Huang's 
thesis, there are nine unknowns; three poles and six zeros. 
The conf iguirat ion of the compensators with these unknowns 
is shown in Figure III-l, v/here the p^ are the poles and 

the z. are the zeros. 

3 

2 . The Plant Equations : 

The plant equations in Figure II-2 as excerpted from 
Reference 4 are almost in state variable form already except 
for the Ci term in equation (3) . By substituting equation 
(2) into equation (3) and collecting terms the follov;ing 
state equations result: 
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u = 



X u 
u 



+ X^a 



- g e 



Z u 

U 



Z 0^ 



a = 



1-2 



1 + Z 



V ' 
o 



1 - z. 

a 



^6^’6t 

V (1-Z.) 

o a 



V (1-Z.) V (1-Z.) 
o a o a 



q = 



M. Z 
a u 

M + 

u X(l-2-) 

o a 



u + 



M. Z 

M + , 
a 1-Z . 

a 



ct + 



M. a-^~) 

ct V 



M + ^ 

^ 1-Z. 

a 



M. + 



M.^6t 

g ^ 






+ 



M. Zr 

"6 " vTT^ 

f o a 



6f + 



M 



M. Z. 
a o e 



6 V (1-Z.) 
e o a 



0 = q 

Substitution of the force and moment coefficients 
from Reference 4 into these equations results in the follow- 
ing plant state equations: 

u = -.0670U + 22.2a -32.29 + .06706^ -6.126^ (5) 

d = -.004017U - .7623a + .9479q -.0013556^ -.15276^ -.07526 (C) 

^ t f e 

q = .003785U - 1 . 826a-l . 452q -.0017326^ +.099536^-1.94516^ (7) 

e = q (8) 

3 . The Actuator Equations 

a,- Thrust actuator: The thrust actuator, in the no- 

tation of Figure II-2, has y^, as input and 6^ as output 
with transfer function: 



1.21 

+ 1.98s + 1.21 
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In differential equation form this becomes: 

I 

1.986^ + 1.216^ = 1.21 y, 
t t t ■'1 

If we let 6^ = 6^^ the state equations for the 

thrust actuator become: 



td 



(9) 

6^^ = -1.216^ -1.986^^ f 1.21 (10) 

b. Flaps actuator: The flaps actuator has input 

Y ^ and 6^ as output with transfer function: 

1.0 

s + 1 . 0 

In differential equation form this becomes: 

6f + 6f = y2 

Thus the state equation for the flaps actuator is 

( 11 ) 

c. Elevator actuator: The elevator actuator has 

input y^ and output 6^ with transfer function: 

121.0 

s^ + 6.6s + 121.0 

The differential equation for this transfer func- 
tion is: 

6 + 6.6 6 + 121.0 = 121. Oy, 

e e e 3 

Letting 6^ = 6^^ the state equations are: 



= -«£ + V2 



6 = 
e 



ed 



( 12 ) 
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(13) 



6 



ed 



-6.6 5 ,-121.05 + 121. Oy,, 

ed e j 



4 . The Compensator Equations : 

In the notation which follows, the ij subscripts 
are defined as follows: the i subscript refers to the 

thrust, flaps and elevator compensators as 1,2,3 respectively. 
The j subscript is the index of the coefficients in each 
compensator. The method of transforming the transfer func- 
tions into state equations when the resulting differential 
equation results in derivatives of the output is as set out 
in reference 5. In the notation of that reference, given an 
n^^ order differential equation: 

x^^^(t) + a,x^^"^^(t) + a iX^^\t)--+ a x(t) = 

1 n-1 n 

b^u^”^ (t) + (t) + ^n-1^^^^ 



where x^^^ (t) = 

dt^ 

A set of state variables for this equation can be de- 
fined as: 



+b (t) 
n 



x^(t) = x(t)-B^u(t) 
X 2 (t) = Xj (t) -B^u (t) 



X (t) = X , (t) “B„ 1 u (t) 
n n-1 n-1 
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Where 



B = b 

0 o 

B, = b, -a, B 

1 1 1 o 

^2 ^ 



n n 1 n-i 



-a n Bi -a B 
n-1 1 no 



With this choice of state variable, the state equations are: 



x^(t) = 






+3j^u (t) 



X2(t) = 



X3<t) 



+B2U(t) 



Xn_i(t) = x^(t) +B_^_ 3 U(t) 

'Vn 

a. Thrust compensator: The thrust compensator has 

input e^, output and the transfer function is: 

= 1000(s+z^) (S+Z 2 ) (s+z^) 

^1 s (s+p^) (s+p^) 

Multiplying out and transforming to a differential 
equation we have: 

y + (p^+p^) yj^+p^P2’y = 1000 e*^+(z^+Z2+Z2)ej^+(z^z^ + Z2Z2 + Zj^Z2)e^ 

+ (Zj^Z2Z2)Sj^ 

Let: = p^+p^ , . a^j = 0 

b^Q = 1000 , bj^j^ = 1000 (z^+Z2-^-z^ ) , b^2~1000 ( z^Z2+Z2Z2+^Z2 ^ 

ba -^ = lOOOz z z 
12 3 
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Thus the state variables are; 



Where 







^13 ^13"^12®ll"^^12 


and the 


state equations are: 






^^2 ^11^1 


(14) 


W 2 


" ^3“^ ^2®1 


(15) 


w^ 


^ ~^12'^2"^13'^3‘'’^13®1 


(16) 



The required output , when these equations are solved 

will then be: ~ ^10^1 

b. Flaps compensator; The flaps compensator has 
input 62 / output Y 2 and its transfer function is: 



Y 2 -400 (s+z^ ) 




Multiplying out and transforming to a differential equation 
we have: 

y^ = -400e2"400z^e2 
Let ^21 ~ ^ 

b 2 Q = -400 , b 23 ^ = -400z^ 

Thus for this first order subsystem, the state variable is: 

“4 ' ^2 ■ ®20®2 

Whore = bjji = -400 , 
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and the state equation is: 



^^4 ^21^2 



( 17 ) 



The required output on solution of this equation 

will be: 

^2 '^4 ^ 20^2 

c. Elevator Compensator: 

The input to the elevator compensator is e^ 

and its output is and it has transfer function: 

-20 (s+z^) (s+Zg) 

~ s (s+p^) 

Multiplying out and transforming to a differential equation 
we have : 



'"’s+Pji-j = -20 e3+{Z5+Zj)e3+Z5ZjC3 



®31 = P 3 ' ^32 = ° 

b3o = -20 , b33 = -20(Z5 +Zj) , bj^ = -2OZ5Z5 

The state variables are: 

Where S2Q=b2Q=-20 , ^30^31 ' ^32"^^32“^31^3l“ ^30^32 

and the state equations are 

(18) 

(19) 

On solution of these equations the required output y^ is 
obtained from: 

^3 ' V®30®3 



*5 = '■'6^®31®3 



w 



^ ~^3l'^6'''^32®3 
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5 . The Feedback and Cost Function Equations 



The feedback loops shown in Figure II-2 are included 
by simply adding three more equations rather than substituting 
in the state equations. This allows for more ease and flex- 
ibility in computer programming of the system model. The 
three equations are: 




The general form of the cost function J discussed in Sec- 
tion above as applied to this system is: 

tf 

j = f (u .-u)^ + (a -a)^ + (0 .-0)^ dt 

v* ref ref ref ” 

o 

6 . The Complete System State Equations : 

Changing the variables in equations 5 to 19 inclusive 
to a consistent set of variables as shown in Table III-l the 
following complete system equations result: 

= -. 067Xj^+22 . 2 x 2~32 . 2x^-6 . 12x.^+ . 067Xg 

x„ = -. 004017x--.7263x.,+ . 9479x--.001355Xt,-.1527xT-.0752x„ 
2 1 2 3 5 7 8 

x^ = .003755x^-1.826x2-1. 452x2". 001732x^+. 09953x^-1. 945xg 




*5 = ^6 

X- = -1 . 21xc-l . 98x,+l . 21x, -+1210 . Oe, 
6 5 6 10 1 
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X? = 

Xg = Xg 

Xg = -121. OXg-6. 6 X 5 + 121 . 0 X 3 ^- 2420.003 

*10 ^ ^ll'*’^!!®! 

*11 ^ ^ 12 ‘‘‘^i 2 ®l 

*12 """^ 12 ^ 11 ^^ 11 ^ 12 ‘^^ 13®1 

^13 " ^ 21®2 

*14 ^15'*’^31®3 

*15 "" "^31^15^32®3 

ei = r^-x^ 

®2 = ^ 2"^2 
®3 "" ^3"^4 

These equations and the defining equations for 
the coefficients in terms of the free parameters (the poles 
and zeros) are coded in Fortran and appear in the sample pro- 
gram appended. The equations appear both in the function e 
valuation subfunction to BOXPLX for calculating J and in 
subroutine GRAPH for plotting the results. Further computer 
symbols used in this programme are explained in Tabxe III-l. 
The computational flow chart showing the inter-relationships 
in the total programme is shown in Figure III-2. 
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Variable 


TABLE III - 1 
Variable Symbols 
State Variable 


Computer Symbol 


u 


^1 


X(l) 


u 




XDOT (1) 


a 


^2 


X(2) 


a 


*2 


XDOT (2) 


q 


^3 


X(3) 




^3 


XDOT (3) 


e 


^4 


X(4) 


0 

e 




XDOT ( 4 ) 




^5 


X(5) 




^6 


X(6) 


6 

t 


"6 


XDOT ( 6 ) 


6 

f 


^7 


X(7) 






XDOT (7) 


6 

e 


^8 


X(8) 


p 

, ^ed 


^9 


X(9) 




^10 


X(10) 


"^2 


^11 


X(ll) 


'^3 


^12 


X(12) 




^13 


X(13) 




^'14 


X(14) 


'"e 


^15 


X(15) 


J 




X(16) 


Pi 




PA(1) 


P2 




PA(2) 


P 3 




PA(3) 


^1 




PA (4) 


^2 




PA (5) 


^3 




PA(6) 


^4 




PA(7) 


^5 




PA{8) 
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TABLE III - 1 (cont,) 



Variable 


Computer Symbol 


^6 


PA(9) 




El 




E2 




E3 




R1 


^2 


R2 


^3 


R3 


Upper bounds 


BU(I) 


Lower bounds 


BL(I) 


Starting values 


XS(I) 


Desired system output 


XDATA(I), XDATAJ(I) 


Intermediate functions of parameters 


AIJ, BIJ, BETAIJ 


Time 


T 


Time step 


DT 
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FIGURE II1"2 Computational Flow Chart 
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D. INVESTIGATIVE APPROACH 

Since it was not known how well BOXPLX would perforin in 
this task of finding the poles and zeros it was decided to 
do all the initial investigation with the known original 
response with the poles and zeros as found by Huang. When 
this phase was completed, it was intended to specify an ar- 
bitrary but reasonable output response and attempt a matching 
of system response to this arbitrary one. 

As was noted earlier the only input used for all phases 
was a step input r2^-~0.25 radians corresponding to a change 
in the pitch angle 0=-O.25 radians. Thus there are two out- 
puts, u and a which are not commanded by a corresponding 
input; and the question arises as to which output to make 
J a function of. This would depend on the situation in 
the real world. If the pitch angle were important due to a 
requirement to stay on a glide slope, then J should be a 
function of 6 . If speed were important due to the distance 
to touchdown, then J should be a function of u . Ideally 
J should be a function of all three outputs, but it may 
not b^ easy to specify just what they should be. The ap- 
proach used in the initial phase was to let J be a function 
of an output not commanded by an input (i.e., speed), then 
let J be a function of pitch angle and finally a function 
of all three outputs. It should be remembered that in this 
initial phase, all three responses, and the poles and zeros 
which produced them, are known. The only other question to 
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be resolved is that of the starting values and upper and 
lower bounds for the variables. Until the capabilities 
of BOXPLX had been checked it was decided to offset the 
starting values approximately 25% from the known values 
(the so-called standard start point) with arbitrary close 
upper and lower bounds. These are shown in Table III - 2 
together with the actual values. 





Pole 1 


Pole 2 


Pole 3 








Actual 


4.0 


10.0 


10.0 








XS 


3.5 


11.0 


9.0 








BU 


5.0 


12.0 


12.00 








BL 


3.0 


8.0 


8.00 










Zero 1 


Zero 2 


Zero 3 


Zero 4 


Zero 5 


Zero 6 


Actual 


1.0 


0.5 


0.5 


0.5 


0.3 


0.4 


XS 


0.75 


0.55 


0.35 


0.55 


0.25 


0.45 


BU 


1.5 


0.7 


0.7 


0.7 


0.4 


0.5 


BL 


0.5 


0.3 


0.3 


0.3 


0.2 


0.3 



Standard Starting Values 
TABLE III - 2 

If BOXPLX proved capable of locating the correct poles 
and zeros to a desired degree of precision, the next step 
was to have it try with the so-called arbitrary start point. 

For this case it will be assumed that the pole and zero lo- 
cations are unknown and the starting values will all be set 
at XS = 10.0, the upper bounds at BU = 20.0 and the lower 
bound at BL = 0 . 0 . 

Assuming success in this initial phase it was planned to 
then specify an arbitrary pitch angle response using Whiteley's 
standard form of response to a step input v;il.h zero position 
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error and 10% maximum overshoot with an appropriate 
(Ref. 7). The speed response would then be checked and if 
"satisfactory", success will have been achieved. If not, 
a speed response using an overdamped second order impulse 
response would be combined with the pitch angle response to 
make a combined cost function. 
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V. RESULTS AND COMMENTS 



A. INITIAL PHASE 

The initial phase, it will be recalled, consisted of in- 
vestigating the capabilities of BOXPLX when the poles and 
zeros of the desired response were known, i.e., the original 
poles and zeros as determined by Huang. It was also during 
this phase that various methods v;ere used to reduce the 
amount of computer time involved in a run and to get an 
idea of how large the value of the cost function could be 
and still have an acceptably close response match. 

Table IV - 1 shows the results of various runs made with 
J a function of speed (i.e., a non-commanded output) and 
starting with the standard start point. The value of the 
cost function varies from J = .00036 to J = .052 . The 
response curves for the smallest value are so superimposed 
as to be indistinguishable and those for the largest value 
are shown in Figures IV - 1, IV - 2 and IV - 3 for u , and 
a and 6 respectively. Although the largest cost function 
is two orders of magnitude larger than the smallest one, the 
response matching is quite acceptable. Thus it appears that 
it is not necessary for BOXPLX to reach a minimum of the cost 
function to give good results and recognition of this fact 
could cut the computer time dramatically. 

The big differences in Table IV - 1 are in the computer 
time involved with each run. It was during this initial phase 
that it was decided to halve the simulation time to 12.5 
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seconds for the calculations in subfunction FE(X) and the 
results show that the computer time is also halved. Another 
ploy attempted was to have BOXPLX set up each complex with 
only n vertices instead of 2n vertices and this again 
halved computer time. However, this last method was found 
to work only when the starting values of the variables were 
close to the minimum and thus would be used only in the 
latter stages of an investigation. 

The next step was to investigate the effect of using the 
standard start point with J still a function of speed but 
enlarging the feasible region by setting the upper bounds at 
20.0 and the lower bounds at 0.0. BOXPLX found a cost func- 
tion of 0.05995 in 600 trials which is about tv/ice the num- 
ber of trials it took to find the last entry in Table IV - 1 
which has a cost function approximately the same. The res- 
ponses with poles and zeros found on this run are shown in 
Figures IV - 4 , IV - 5 and IV - 6 , for speed angle of attack 
and pitch angle respectively. • Again the results are quite 
acceptable. Another interesting result of this run is that 
the pole and zero values differ, in some cases appreciably 
both from those obtained from the last run in Table IV - 1 
and from the actual values, and yet the responses are quite 
similar. The values obtained from this run are shown belov; 
in Table IV - 2 . A further mention of this sensitivity of 



the results will be made later. 
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Pole 1 



Pole 2 



Pole 3 



Found 


4.409 


0.437 


10.013 








Actual 


4.0 


10.0 


10.0 










Zero 1 


Zero 2 


Zero 3 


Zero 4 


Zero 5 


Zero 6 


Found 


0.664 


0.336 


1.107 


0.447 


0.237 


0.443 


Actual 


1.0 


0.5 


0.5 


0.5 


0.3 


0.4 



Parameters found, standard start, wide bounds, J=f (u) = .05995 

Table IV - 2 

Next, the cost function J was made a function of the 
pitch angle output, 0 (the commanded ouput) again using the 
standard start point and the larger feasible region. In 600 
trials, taking 17 minutes BOXPLX found a J=6.8 x 10 ^ and 
the results of the simulation of the poles and zeros found 
are shown in Figures IV - 7 , IV - 8 and IV - 9 for u , a and 
9 respectively. Although the cost function as a function 
of 0 cannot be compared directly v/ith the one which is a 
function of speed because of the scale factors involved it 
can be considered to a first approximation to be approximately 
three orders of magnitude smaller. Thus, this J as a func- 
tion of 6 would be roughly equivalent to one as a function 
of u of .068. But the important point to note is that the 
response match is much better than when J was a function of 
u and in fact the values found for the poles and zeros cor- 
respond to the actual ones with a maximum of 12% with another 
8.8% and the rest less than 2%. Thus it would appear that 
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better convergence is obtained when J is made a function 
of the commanded output. 

Next came the real test of the capabilities of BOXPLX. 
Here it was assumed that the actual poles and zeros of the 
response are unknown and BOXPLX would start from a so-called 
arbitrary start point in a large feasible region. For this 
purpose, as mentioned earlier, the upper bounds were all set 
at 20.0, the lower bounds at 0.0 and the starting values at 
10.0, as might be assumed by a reasonably competent engineer, 
examining this model. The cost function was made a function 
of speed and then a function of the pitch angle, as before, 
to check on the finding in the preceeding paragraph. 



First, with 


J = f(u) 


, BOXPLX 


found J 


= 0.99016 


in 736 


trials 


with the 


responses 


shown in 


Figures 


IV - 10, 


IV - 11 


and IV 


- 12 for 


u , a and 6 respectively 


and found 


the poles 


and zeros shown 


in Table 


IV - 3 below. 








Pole 1 


Pole 2 


Pole 3 








Found 


5.328 


8.411 


12.289 








Actual 


4.0 


10.0 


10.0 










Zero 1 


Zero 2 


Zero 3 


Zero 4 


Zero 5 


Zero 6 


Found 


0.229 


0.722 


1.342 


1.317 


0.084 


1.130 


Actual 


1.0 


0.5 


0.5 


0.5 


0.3 


0.4 



Parameters found, arbitrary start, J = f (u) = 0.99016 

Table IV - 3 
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It will be noted that even though these values differ 
markedly from the actual ones, the u and ot responses are 
good. But, there is a slowly decreasing steady state error 
in the response due to the very small value of Zero 5 which 
has almost cancelled the pole at the origin in the pitch an- 
gle compensator. 

In the hope of improving these results, the programme 
was rerun with the values in Table IV - 3 as the starting 
point, with the same upper and lower bounds (the bounds could 
be changed to narrow the feasible region, but it was decided 
to see how well the convergence went without doing this) . 

In a further 555 trials or a total of 1291 BOXPLX came up 
with a J = .08477 and the values shown in Table IV - 4 below. 





Pole 1 


Pole 2 


Pole 3 








Found 


8.628 


5.341 


11.150 








Actual 


4.0 


10.0 


10.0 










Zero 1 


Zero 2 


Zero 3 


Zero 4 


Zero 5 


Zero 6 


Found 


0.3000 


1.589 


0.497 


0.870 


0.000 


0.912 


Actual 


1.0 


0.5 


0.5 


0.5 


0.3 


0.4 


Parameters found. 


arbitrary start. 


J = f (u) 


= 0.0848 










Table IV 


- 4 







The responses of the system with these poles and zeros 
are shown in Figures IV - 13, IV - 14, and IV - 15 for u , a 
and 0 respectively. 
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It is to be noted that, although J is smaller by an or- 
der of magnitude and the u and a responses are well nigh a 
perfect match, there is now a steady state error in the 9 
response due to a cancellation of the pole at the origin in 
the pitch angle compensator. This situation is not likely 
to improve with further trials with J a function of speed. 
This is the first indication of a worse response with more 
trials which will arise again later. 

The same arbitrary start point was used on the next run 
with J a function of 0 (the commanded output) . This 
time, the results were more encouraging. In 900 trials with 
J - 4.1 X 10 ^ (again not directly comparable with the value 
when J is a function of speed) BOXPLX came up with the 



values 


shown it 


Table IV 


- 5. 










Pole 1 


Pole 2 


Pole 3 








Found 


5.259 


8.708 


9.981 








Actual 


4.0 


10.0 


10.0 










Zero 1 


Zero 2 


Zero 3 


Zero 4 


Zero 5 


Zero 6 


Found 


0.467 


0.249 


1.682 


0.447 


0.325 


0.374 


Actual 


1.0 


0.5 


0.5 


0.5 


0.3 


0.4 


Parameters found 


, arbitrary start. 


J = f( 


) = 4 . 1 X 


10-6 








Table IV 


- 5 






The 


responses for u 


, a and 


0 are 


plotted in 


Figures 


IV - 16 


, IV - 17 


and IV 


- 18 respectively. 


The match 


is near 


perfect 


for a 


and 0 


and quite 


acceptable for u . 


This 



run reinforced the earlier assertion that better results 
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are obtained when J is a function of the commanded output 
(at least for this system) . 

One new wrinkle which appeared during this last run was 
the case where BOXPLX was not able to find a feasible ver- 
tex at a given trial and restarted to compute at the best 
vertex. However, during the course of computation after 
this restart, it was again not able to find a feasible ver- 
tex and the best vertex was still the one at which it had 
restarted. Thus exit from the subroutine occured. The re- 
medy in this case was to change the value of R , the first 
random number of the sequence generated in developing the 
initial complex of vertices. It must be admitted that this 
is probably a matter of luck as it is with all random numbers, 
but the results of the run speak for themselves. 

The last test in this initial phase was to make J a 
function of all three outputs as an example of this admit- 
tedly artificial situation, since it is not likely that, for 
a given problem, one would know what all the outputs should 
look like, much less specify them. 

The problem arising here is to put weighting factors in 
the cost function so that the response variations from the 
desired response may be penalized depending on the importance 
of one or more outputs over others. If all are considered 
equally important which would normally be the case if they 
were specified it remains only to balance the cost function 
component variations due to scale factors. This was fairly 



64 



I 

I 

I 

1 



’ .^1 

- ~ 4 



- s 



J 



: *> ? 



r. 



: i 



3 

% 




I 



easy in this case because of the data from past runs. The 

value of J for the arbitrary start point when J was a 

function of speed was 0.9989046 x 10^^ and the equivalent 

47 

for pitch angle was 0.7685555 x 10 .A trial run with J 
a function of all three outputs without weighting factors 
produced J = 0.9990991 x 10^^. Subtracting the first two 

48 

from the last resulted in a J for alpha = 0.1176444 x 10 
Thus as a rough approximation, the speed component of the cost 
function was made 0.001, the angle of attack component, 0.1 
and the pitch angle component, 1.0. 

With this cost function, the best vertex after 1120 trials 
with J = .00042 was as shown in Table IV - 6 . 





Pole 1 


Pole 2 


Pole 3 






Found 


6.339 


7.203 


10.597 






Actual 


4 . 0 


10.0 


10.0 








Zero 1 


Zero 2 


Zero 3 


Zero 4 


Zero 5 Zero 6 


Found 


1.003 


0.252 


1.098 


0.641 


0.240 0.524 


Actual 


1.0 


0.5 


0.5 


0.5 


0.3 0.4 


Parameters found. 


arbitrary start. 


J = f (u ,a 


,0) = 0.00042 








Table IV 


- 6 




The 


system responses 


for these 


poles and 


zeros are plot- 



ted in Figures IV - 19, IV - 20 and IV - 21 for u, a and 6 
respectively. The response matching is almost perfect al- 
though the pole and zero values differ considerably from the 
actual values. 
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B. ARBITRARY DESIGN PHASE 

At this point, it was considered that the intricacies 
of BOXPLX and the factors affecting desired results were 
understood well enough to form a strategy by which a multi- 
variable system of a given configuration could be designed 
to have any time response desired to any input. 

The strategy was formulated as follows: 

a. Specify the time response and make the cost func- 
tion a function of the output which is commanded by the in- 
put. 

b. Let BOXPLX try to match that output response for 
600 to 900 trials and examine plots of all outputs. 

c. If other outputs are reasonable, success has been 
achieved. If not, specify another output which may be import- 
ant and make J a function of both outputs. This other out- 
put specification should be a compromise between what the 
system seems capable of achieving and a reasonable desired 
response . 

d. Continue iteration until all responses are satis- 
factory. 

Consequently, for this system it was decided to specify 
the pitch angle response to a step input of -0.25 for r^ as 
shown in Figure IV - 22. This particular step response was 
obtained from a state variable simulation of Whiteley's 
standard form for a step response with a maximum of 10% 
overshoot on zero position error (Ref. 7) with an = 8.75. 



69 



pitch angle - radians 



+ .05 



0.0 



-0.05 



- 0.10 



-0.15 



- 0.20 



-0.25 



-0.30 

























1 








1 








1 

i 










■ 1 

1 












- - 








i 


1 











10 15 20 



t iir.e-soconds 
Figure IV - 22 



70 








m 






r*4 



• 0 



I 

( 

a*- 



I 






The value of w was chosen from an examination of the nat- 

o 

ural time constant of the system as revealed by previous data, 
which showed a rise time of approximately 0.5 seconds. Ini- 
tials runs with = 2.5 which gave a rise time of 1.75 se- 

conds showed that BOXPLX was having a difficult time trying 
to design the compensator for such a slow response. The 
values for Pole 3 kept ending up on the upper bounds, and 
they were raised twice to 30.0 and 40.0 before an = 8.75 

was chosen. The computer program at the end of this paper 
shows how the response was simulated and cards punched to 
act as reference input to the function minimization program. 
This response is shown in Figure IV - 22. 

Thus a run was made with this response a reference and 
J function of 0 . In 900 trials with J = . 00156/ BOXPLX 
found the poles and zeros shown in Table IV - 7 . 



Pole 1 


Pole 2 


Pole 3 








3.669 


14 . 644 


19.890 








Zero 1 


Zero 2 


Zero 3 


Zero 4 


Zero 5 


Zero 6 


0.802 


0.397 


0.034 


2.054 


0.752 


1.220 



Parameters found, arbitrary start, J = f(6) = .00156, ar- 
bitrary response 

Table IV - 7 

The response curves of the system with the above poles 
and zeros are plotted in Figures IV - 23, IV - 24 and IV - 25 
for u , a and 0 respectively. The theta response matches 
the reference arbitrary response fairly well and the alpha 
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response looks reasonable but there is an extremely long 
settling time for the speed to return to zero. This was 
due to the extremely small value of Zero 3 in the speed com- 
pensator which almost cancelled its pole at the origin. 

In the hopes that more trials might improve the results 
the programme was run again starting from the values in Ta- 
ble IV - 7. After a further 900 trials or 1800 tota], with 
J = .00088, BOXPLEX found the values give in Table IV - 8 
below: 



Pole 1 


Pole 2 


Pole 3 








7.173 


8.635 


14.933 








Zero 1 


Zero 2 


Zero 3 


Zero 4 


Zero 5 


Zero 6 


0.884 


0.243 


0.000 


0.000 


0.809 


0.810 



Parameters found, arbitrary start, J = f(0) = .00088, ar- 
bitrary response. 

Table IV - 8 

The results of a system simulation with these values are 
shown in Figures IV - 26, IV - 27 and IV - 28 for speed, an- 
gle of attack and pitch angle respectively. It can be noted 
that the match to the arbitrary theta response is excellent 
and that the short period oscillations both in it and in the 
alpha response have been reduced in amplitude and duration. 
But, the speed response now has a constant steady state error 
due to the fact that decoupling has been lost. It will be 
recalled that decoupling of this multi-variable system re- 
quired a pole at the origin in each of the compensators, but. 
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since Zero 3 in the speed compensator went to the lower 
bound of 0.0, it effectively cancelled out this pole. 

Consequently, it was decided to change the lower bounds 
for all the zeros to an arbitrary value of 0.25 and rerun 
the programme from the arbitary start point. This time, 
BOXPLX found the following values in Table IV - 9 after 
900 trials with J = .00180. 



Pole 1 


Pole 2 


Pole 3 








8.994 


19.519 


14.510 








Zero 1 


Zero 2 


Zero 3 


Zero 4 


Zero 5 


Zero 6 


0.331 


0.298 


2.332 


0.763 


0.942 


0.702 



Parameters found, arbitrary start, J = f(6) = .0018, arbi- 
trary response, BL = .25 on all zeros. 

Table IV - 9 

The response of the system with these values for the 
poles and zeros are plotted in Figures IV - 29, IV - 30, and 
IV - 31 for u , ct and 0 respectively. The results show 
some secondary overshoot on the 0 response, some increased 
oscillations in the a response but the speed response is 
beginning to look like what one would expect from the con- 
ditions . 

To check if further trials would improve the results the 
programme was rerun for a further 900 trials starting from 
the values in Table IV - 9 . The results of these 1800 trials 
with J = .00099 are shown in Table IV - 10 below. 
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Pole 1 



Pole 2 



Pole 3 



7.007 8.231 15.372 

Zero 1 Zero 2 Zero 3 Zero 4 Zero 5 Zero 6 

0.240 0.253 0.567 0.240 0.842 0.828 

Parameters found arbitrary start, J = f(0) = .00099, arbi- 
trary response. BL = 0.25 on all zeros. 

Table IV - 10 

The responses of the system with these poles and zeros 
are shown in Figures IV - 32, IV - 33 and IV - 34 for u , 
a and 0 respectively. As was mentioned before as being 
sometimes the case, these further trials did not improve the 
theta and alpha responses significantly but definitely wor- 
sened the speed response to one with a long settling time 
and some long period oscillation. 

Thus it was now time to backtrack to the previous run 
and specify a suitable speed response and make the cost func- 
tion a function of both u and 0 . The response chosen in 
concert with colleagues cognizant in this area was the one 
shown in Figure IV - 35. This response was obtained from a 
state variable simulation of the impulse response of a second 
order overdamped system as shown in the computer programme at 
the end of this paper. 

After revising the function minimization programme to the 
final form shown at the end of this paper as a sample of all 
such programmes, to include a balanced J as a function of 
u and 0 , BOXPLX , in 1313 trials found the pole and zero 
values given in Table IV - 11. 
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Pole 1 



Pole 2 



Pole 3 



12.236 
Zero 1 



22.236 
Zero 2 



15.354 
Zero 3 



Zero 4 



Zero 5 



Zero 6 



1.444 0.302 1.305 0.360 0.436 1.290 

Parameters found, arbitrary start, J = f(u,0) = .04667, ar- 
bitrary response BL = 0.25 on all zeros. 



Table IV - 11 



The results are shown plotted in Figure IV - 36, IV - 37 
and IV - 38 for u , a and 0 respectively and the excellent 
results can be noted. This was also the best that BOXPLX 
could do since it exited after two restarts with no improve- 
ment in the best vertex. One need only glance at the ori- 
ginal response of the system to note the improvement obtained. 



C. SENSITIVITY OF RESULTS 

As was noted earlier, some excellent response matching 
has occured in the process of this investigation where, in 
the initial phase, the poles and zeros found by BOXPLX varied 
widely from those which produced the reference response. As 
a result there could be orders of magnitude difference in 
values of J and still the results could be acceptable. 

Thus, of itself, a "small" value for J is not necessarily 
indicative of a good response matching. 

A case in point is that where J was a function of all 
three outputs. The graphed responses were so superimposed 
as to be almost indistinguishable yet the value of Pole 1 
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varied from its actual value by over 50% and that for Zero 3 
by over 100%. Both of these parameters were in the speed 
compensator and it might be said that these variations to- 
gether with those of the other three parameters in the same 
compensator might have interacted with each other to produce 
an overall acceptable response. Nonetheless, for all the 
parameters, the results show insensitivity for variations 
in the parameters of up to 25 and 30%. 
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V. CONCLUSIONS & RECOMl*lENDATIONS FOR FURTHER RESEARCH 



A. CONCLUSIONS 

The investigations carried out in this thesis result 
in the following conclusions: 

a. It is possible in a multivariable, steady-state 
decoupled system with cascade compensators and unity feed- 
back to design the compensators by computer methods to have 
the system match any reasonable specified response to a 
given input. 

b. Where computer time is a problem, there exist 
several methods, as outlined in this paper, for the function 
minimization subroutine used here to reduce the time taken 
to produce an acceptable response match. 

c. Where only a limited number of parameters are 
available, any arbitrary desired response should have time 
constants compatible with the natural time constants of the 
system. 

d. Where possible, a cost function which is a func- 
tion of two or more outputs, when balanced, will produce 
better results than when it is a function of only one out- 
put. 

e. Where the cost function is a function of only 
one output, the output to be chosen should be the one which 
is commanded by the input. 

f. The results show insensitivity of response to 
variations in parameters of up to 25 and 30%. 
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B. RECOMMENDATIONS FOR FURTHER RESEARCH 



The following areas are recommended to extend the use 
of the design tool developed in this paper: 

a. The gains in the compensators, which were kept 
constant in this investigation, could be made free parameters 
to extend the flexibility of the design procedure. 

b. The numbers of poles and zeros in each compen- 
sator were left at those developed in Huang's paper. An in- 
vestigation should be made to see if the response matching 
could be improved by increasing the numbers of poles and 
zeros, and by how many. 

c. The flexibility of the method would be further 

increased if provision were made to include complex poles 

and zeros as well as real ones. By using transfer functions 
s^-f a + b 

of the form where a and b are real num- 

(s+p^) (S+P2) 

bers, the parameters would be a combination of poles, zeros 
and coefficients. 

d. Military specifications for aircraft contain 
specific values for limits on parameters of aircraft motion 
such as short period and phugoid frequency and damping. If 
the equations were altered to obtain these parameters of the 
aircraft response in terms of the compensator poles and zeros, 
a matching to military specifications could be obtained by 
insertion of a suitable implicit constraint evaluation sub- 
routine to BOXPLX. 
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APPENDIX A 



DESCRIPTION OF SUBROUTINE BOXPLX 

BOXPLX is a Fortran-63 subroutine which finds the mini- 
mum of a general non-linear function within a feasible region 
by the complex method of M.J. Box. Explicit constraints are 
defined as upper and lower bounds on the independent variables. 
Implicit constraints may be arbitrary functions of the vari- 
ables. Two function subprograms to evaluate the cost function 
and the implicit constraints respectively msut be appended 
by the user. 

The complex method is an extension and adaptation of 
the simplex method of linear programming. Starting with any 
one feasible point (XS(I)) inn-dimensional space for n vari- 
ables, a "complex" of 2n vertices (sets of variables) is 
constructed by selecting random points (via a random number 
generator) within the feasible region. For this purpose, 
n coordinates are first randomly chosen within the space 
bounded by the explicit constraints. This defines a trial 
initial vertex, which is then checked for possible violation 
of the implicit constraints. If one or more are violated, 
the trial initial vertex is displaced half of its distance 
from the centroid of previously selected initial vertices. 

This displacement process continues for up to 5n displace- 
ments until the vertex has become feasible. If not, exit 
occurs from the subroutine. 
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If an initial complex is established, the basic compu- 
tational loop is initiated. The current worst vertex, i.e., 
the one with the largest value of the cost function, is re- 
placed by its over-reflection through the centroid of all 
other vertices. (If the verte.x to be replaced is considered 
to be a vector in n-space, its over-reflection is opposite 
in direction, increased in length by a factor 1.3 and col- 
linear with the replaced vertex and the centroid of all other 
vertices) . 

When this overreflection is not feasible or remains the 
worst vertex, it is displaced half-way toward the centroid. 

If NV displacements occur without a successful result, it is 
over reflected through the best vertex, i.e., the one with 
the smallest value of the cost function. If a successful 
result is still not achieved, the whole process restarts at 
the best vertex or the centroid, whichever has the smaller 
cost function. 

However, in most cases, all vertices converge to a point 
which is taken as a solution. Convergence is considered at- 
tained when 2n consecutive cost function values have been 
within 10 of a base value. 

The parameter list for BOXPLX is flexible and allows the 
user to choose: 

a. the number of variables 

b. the number of auxiliary variables (implicit con- 
straints) 

c. the frequency of printed output 

d. the number of trials allov/ed 

e. the first random number to be used 
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f. the starting values 

g. the upper and lower bounds 

If BOXPLX does not find a minimum cost function at the end 
of the number of trials specified it prints out the values 
for the best vertex reached up to that point. If further 
trials are considered necessary, these values should be used 
as the starting values for the next run. 

Obviously, the starting values selected must lie in the 
feasible region, i.e., between the upper and lower bounds, 
or BOXPLX. will not even start. 

If calculation restarts at the best vertex after being 
unable to find a feasible vertex and the same situation oc- 
curs again, the "new" best vertex is compared with the pre- 
vious one. If it is better, calculations restart again at 
this "new" best vertex, otherwise exit from BOXPLX occurs 
with the previous best vertex being printed out. 
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THE SIMULATION ANO PLOTTING SUBROUTINE 
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PLOTS THE RESPONSE OF THE CONPENSATED SYSTEM AND THE DESIRED RESPONSE 
SUBROUT I NE GR APH ( PA, XDATAl ,XDATA2) 

DIMENSION PA(9) ,TIME(501) ,SPEcO(501) ,DS PEED (501) ,XDATA1( 501) , ALPHA 
*(501) , THETA (501) , DTHETA( 501 ) ,XDATA2( 501 ) 

REAL* 8 X( 15) ,XOOT( 15), T,DT, TITLE! 1 2 ) / • M ACN AMAR • ,’A M.A.S.',10*' 
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THE FUNCTION MINIMIZATION SUBROUTINE 
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